BCS-BEC crossover of collective excitations in two-band superfluids 
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We use the functional integral approach to study low energy collective excitations in a continuum 
model of neutral two-band superfluids at T = for all couplings with a separable pairing interaction. 
In the long wavelength and low frequency limit, we recover Leggett's analytical results in weak 
coupling (BCS) for s-wave pairing, and further obtain analytical results in strong coupling (BEC) 
for both two and three dimensional systems. We also analyse numerically the behavior of the out- 
of-phase exciton (finite frequency) mode and the in-phase phonon (Goldstone) mode from weak 
to strong coupling limits, including the crossover region. In principle, the evolution of Goldstone 
and finite frequency modes from weak to strong coupling may be accessible experimentally in the 
superfluid phase of neutral Fermi atomic gases, and could serve as a test of the validity of the 
theoretical analysis and approximations proposed here. 

PACS: 74.40.+k,74.20.Fg,05.30.Fk,03.75.Kk 



I. INTRODUCTION 

A two-band model for superfluidity was introduced by 
Suhl et ali in 1959 soon after the BCS theory in or- 
der to allow for the possibility of multiple band crossing 
through the Fermi surface. Suhl et ali observed that a 
larger number of energy bands crossing the Fermi sur- 
face could increase the overall electron state density and 
lead to the onset of additional interactions. Since then, 
this model has been used to describe high temperature 
superconductivity in copper oxides and more recently it 
has been used in connection to MgBj^. 

In the seminal work by Leggetl4, the existence of col- 
lective phase modes in two-band superfluids has been 
predicted in both neutral and charged systems. Within 
the weak coupling (BCS) limit, Leggett showed that if 
s-wave interactions are attractive in both bands, an un- 
damped long-wavelength exciton (finite frequency) mode, 
as well as an undamped long-wavelength phonon (Gold- 
stone) mode may exist. These modes were further stud- 
ied theoreticallji^& in the BCS limit, however, there was 
no experimental evidence of their existence until man- 
ifestations of two-gap behaviouri2i2ii2iii were found in 
MgB 2 . 

However, it has not been possible to study the evolu- 
tion of collective spectra from weak (BCS) to strong cou- 
pling (BEC) limit until very recently. Advances in exper- 
iments with neutral Fermi gases enabled the tunning and 
control of two-particle interactions between atoms in dif- 
ferent hyperfine states by using Feshbach resonancesi 2 * 1 ^. 
This kind of control is not fully present in standard 
fermionic condensed matter systems, and has hindered 
the development of experiments that could probe sys- 
tematically the effects of strong correlations as a function 
of coupling or density of fermions. It was thought the- 
oretically for many years that a weakly coupled (BCS) 
superfluid could evolve smoothly into the limit of tightly 
bound pairs which undergo Bose-Einstein condensation 
(BEC) 14 i 15 i 16 i 17 i 18 i 19 i 20 . It was not until recently that 
the first experimental evidence that hyperfine states of 
40 K (Ref^iiS 2 .) and 6 Li (Refi 23 i 24 i 25 i 26 ) can form weakly 



and tightly bound atom pairs (Cooper pairs), when the 
magnetic field is swept through an s-wave Feshbach res- 
onance. 

Considering these recent findings in condensed mat- 
ter systems and advances in atomic physics experiments, 
we expand Leggett's calculation of collective modes in 
neutral two-band (s-wave) superfluids to all couplings by 
following a similar one-band approach 27 . These collec- 
tive modes for two-band s-wave systems are undamped 
in the low- frequency and low- momentum limits, provided 
that the two-quasiparticle threshold is not reached. We 
present results of the evolution of the finite frequency 
and Goldstone's modes from weak coupling (BCS) and 
to strong coupling (BEC) limits, and discuss briefly the 
possibility of observing these modes in experiments in- 
volving multi-component ultracold atomic Fermi gases. 
While it is still a matter of debate that extentions of Ea- 
gles—, Leggett'siii, and Nozieres and Schmitt-Rink's 16 
(NSR) suggestions are good quantitative description of 
the crossover phenomena, the evolution of collective 
modes can serve as a test of these ideas. 

The rest of the paper is organized as follows. In sec- 
tion [HJ we discuss the effective action and the saddle 
point approximation for a two-band continuum Hamil- 
tonian with attractive interactions in the s-wave chan- 
nel and a Josephson interband coupling term. We anal- 
yse the effects of Gaussian fluctuations in section IIIII 
and derive effective amplitude and phase actions, from 
which the in-phase (Goldstone) and out-of-phase (finite 
frequency) collective mode spectra are calculated in sec- 
tion IIVI The evolution of these modes is analysed as 
a function of interaction strengths from the BCS to the 
BEC limit both analytically and numerical in section Ivl 
In section lv*Tl we summarize our conclusions and propose 
that ultra-cold Fermi atoms can be used to test our re- 
sults regarding the evolution of the finite frequency and 
Goldstone modes in two-band superfluids. Finally, in Ap- 
pendix A, we present details of the matrix elements in- 
volved in the phase and amplitude effective actions, while 
in Appendix B, we show the long wavelength expansion 
coefficients needed to evaluate the phase modes. 
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II. EFFECTIVE ACTION METHOD 

A Hamiltonian for multi-band (or multi-component) 
superfluids with singlet pairing can be written as 

H = Y ^(k^Lk.cr ™*^ + ^ ^ U n7n (Q)pn,qPm,— q 
n,k,cr n,m,q 

+ VnmiKk'^l^b^^ (1) 

n,m, k,k'.q 

where the indices n and m label different bands (or com- 
ponents), a labels the spins (or pseudo-spins), and k la- 

a n k t i s ^ ne fermion 



In addition, 



bels the momentum 
creation operator, p„, q = a ri,k-q,<T a ™,k, CT is the den- 
sity operator, and 6^ k q = aJ l , k+q/ 2,t Q v l ,-k+q/2,.l is thc 
pair creation operator. Here, £„(k) = e„(k) — /x, where 



e„(k) = £„ :0 + k 2 /2m„ 



(2) 



(with 7i = 1) is the kinetic energy of fermions. The refer- 
ence energies are Si t o = and £2,0 — E > (see Fig.^l; 
m n is the effective mass for the n th band. Furthermore, 
the terms U nm and V nm correspond to the Coulomb and 
pairing interaction matrix elements, respectively. Since 
in this paper we arc interested only in thc neutral case, we 
will ignore Coulomb interactions, and consider only the 
pairing term. This choice is more appropriate to ultra- 
cold atomic Fermi gases, while the inclusion of Coulomb 
terms is more appropriate in the case of superconduc- 
tors. The discussion of the charged case is postponed to 
a future manuscript. For the neutral case, we take 



vw(k,k') = i4mr;(k)r m (k') 



(3) 



to be separable, and we consider in general a two-band 
system with distinct intraband interactions Vn and V22, 
and interband interactions V12 and V21. Notice that the 
off-diagonal terms V\2 and V21 play the role of Josephson 
coupling terms. For the purpose of this paper we will 



J 



consider all V nm to be negative. The T„(k) coefficients 
are symmetry factors characterizing the chosen angular 
momentum channel. 




FIG. 1: Schematic figure of two-bands with reference energies 
£1,0 = for the band 1 and £2,0 = Eo for band 2. 



In the imaginary-time functional integration formalism 
(f'3 = 1/T, h = kg = 1), the partition function is written 
as 



Z = J D[a\a]e' 
with an action given by 



S= dT 
Jo 



E <k,>)(drK, k , ff (T)+ff(T) 
n,k,er 



(4) 



(5) 



The Hamiltonian from Eq. can be re-written in the 
form 



H(r) = £ U^W nM Jr)a nM , a (r) 

n,m,k,cr 

+ Y B, t i (q,T)K m B m (q,r), (6) 



n,m,q 



with B„(q, t) = ^ k r„(k)& nj k jq . We first introduce 
the Nambu spinor ip^ip) = ( a^ n p ^,a n ^ p [ ), where p = 
(k, iwi) is used to denote both momentum and Matsub- 
ara frequencies (wi = (2£+ (3). Furthermore, we use 
thc Hubbard-Stratanovich transformation 



exp{- Y, B n (q)V nm B m (q)} = f D[*\ $] exp { Y ®U<l)9n m $ m (q) + Y [ B n(<z)$«(g) + h.c] } 

n,m,q n,m,q n,q 

r 



to decouple the fermionic degrees of freedom at the 
expense of introducing the bosonic fields $„((?), with 
q = (q, ivi) where vi — 2£ir//3. The tensor g mn associ- 
ated with the bosonic pairing fields $ m (q) can be written 

as 



5n 5i2 
.921 .922 



1 



V22 —V12 
-V21 V n 



(8) 



detV 

where detV = V11V22 — V12V21. Performing an integra- 



tion over the fermionic part (D[i/j^ , ip}) leads to 

S = -0 Y $ n(9).9nm$ m (9) 

n,m,q 

+ Y [Uk)<W -TrlnG- 1 ] . 



(7) 



(9) 



3 



Here, the inverse Nambu matrix is 

g; 1 = $„(- g )r;(^)a_ + $t( g )r„(^)a + 

+ [iwta - £ n (k)o- 3 ] 5 PiP i, (10) 

where a± — (<7i ± <72)/2 and e>"i are the Pauli spin matri- 
ces. 

We choose an approximation scheme where the field 
is written as a sum of a r-independent (stationary) 
part and a r-dependent contribution 

- P 1 ) = $n(q) = A n , <5g,0 + A„(g). (11) 

We first write the inverse Nambu matrix in terms of two 
matrices: 



G n — G n o (1 + G„ o 

G»,l )• 



(12) 



The first one is the saddle point inverse Nambu matrix 
given by 

r -l ( *w ~ £n(k) A; j0 r n (p) \ 
n0 - dpy I A n , r;( P ) lw + e„(k) I > ^ J 



and the second one is the fluctuation matrix 



G n i 1 — 



A*(-q)r n ( 



p+p' - 



A n (q)r*(E±P_) 







(14) 



used in the expansion of the natural logarithm of 
TrlnG" 1 = TrlnG n . - ]T ^TrfGn.oG^r 1 ) 1 - 

i=l 

, (15) 

Within this approximation, we expand the action to 
second order in the fluctuation field A n (q). This proce- 
dure produces an effective action of the form S = S0 + S2, 
where 

Sq = — /3 <jVi, m A* >0 A m ,o 

n,m 
n,p,p f 

is the effective saddle point action, and 

S-2 = Y1 G n , (12)G n , 1 - 1 (23)G n:0 (34)G n , 1 - 1 (41) 

is the Gaussian correction to it. The notation (ij) in 
the G matrices is understood as the momentum labels 
(pi,Pj). An important comment about So is in order. 
Writing A nj o in terms of its amplitude and phase 



A„ j0 = |A„,o| exp(itp n ), 



(16) 



the first term of So becomes the sum of two contributions: 
the standard band-diagonal terms — ^ n g„„| A n0 | 2 , 
and the band-off-diagonal — 2f3gn\ Ai,o| I A2,o| cos(</J2 — 



(fi) corresponding to the Josephson coupling between 
bands. Since g i2 = — Via/detV, with detV > and 
all V nm < 0, the saddle point thermodynamic poten- 
tial f2o = So//3, has its quadratic term of the form 
(|V22||A li0 | 2 + |Vn||A 2 ,o| 2 - 2|Vi2||A li0 ||A 2)0 |cos(¥>2 - 
</3i))/ det V, which shows explicity the Josephson energy. 
For the case chosen, where all V nm are negative, the in- 
phase <p2 — fi is the only stable solution. However, if 
V12 were positive, another stable solution would appear 
where ip 2 — ipi+n. This is the so-called 7r-phase solution, 
which we will not discuss here. 

From the stationary condition dSo/dA^q) = we ob- 
tain the order parameter equations 



n, k' 



tanh ^ g "< k ') 

t r„(kOA», r n (kO 2 * (17) 



where E n (k) = (£ 2 (k) + |A n (k)| 2 ) a is the quasiparticle 
energy spectrum. Note that the order parameter 



A„(k) = A n!0 r„(k) 



(18) 



is a separable function of temperature T and momentum 
k. 

In this manuscript, we focus on s-wave superfluids, and 
thus consider only the zero angular momentum channel 
of the interaction V nm (k, k'). In addition, instead of tak- 
ing r„ (k) = 1 (independent of k) which would cause ul- 
traviolet divergences (logarithmic in two dimensions and 
linear in three dimensions) in the integrations over mo- 
mentum for the order parameter equation, we take 



r„(k) = i/(i + fc/fc«,o) 1/2 



(19) 



as the corresponding symmetry factor in two dimensions 
(2D), and 



r„(k) = 1/(1 + fc/*V.,o) 



(20) 



as the corresponding symmetry factor in three dimen- 
sions (3D). Here k Hi o ~ -R^o, where R n ,o plays the role 
of the interaction range in real space, sets the scale at 
large momenta, and it is necessary to produce the phys- 
ically correct behaviour of a generic interaction at short 
wavelengths. In the case of 2D, V^„(k, k') ~ 1 for small 
k and k'), while V nn (]s., k') ~ 1/ykk' for large k and 
k'), when a generic real potential is used 2 ^. There is 
no ultraviolet divergence in our theory since the mo- 
mentum integrations always produce finite results. This 
choice of interactions has the advantage of making unnec- 
essary the introduction of the T-matrix approximation 
to renormalize the order parameter equation, and rede- 
fine the interaction amplitude in terms of the two-body 
binding energy 29 30 . In the case of 3D, V nn (\s., k') ~ 1 
for small k and k'), while V^ n (k, k') ~ 1/kk' for large 
k and k', when a generic real potential is usecU^. No- 
tice that, our three dimensional interaction T„(k has the 
same behaviour at both low and high momenta with the 
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one used by Nozieres et al. 16 T(k) = 1/^1 + (k/k ) 2 . Ei- 
ther choice (ours or NSR's) produces qualitatively similar 
results as it will be later discussed. 

Furthermore, the order parameter equations need to 
be solved self-consistently with the number equation 



^= E 



i 



gn(k) 

2E n (k) 



tanh 



2k h T 



(21) 



which is obtained from JV = —dflo/dfi, where fio is the 
saddle point thermodynamic potential. The inclusion of 
fluctuations are very important for the number equation 
near the critical temperature of the system, and in this 
case the saddle point thermodynamic potential needs to 
be corrected to fl = Oq + ^G: where tic should be calcu- 
lated at least at the Gaussian level. In this manuscript, 
however, we limit ourselves to fluctuation effects at low 
temperatures, as discussed next. 



block diagonal indicating that the two bands are uncou- 
pled. This corresponds to the case where the Josephson 
coupling V12 = V21 = 0, since 312 = — V\ij det V as indi- 
cated in Eq. JSJ|. 

Performing Matsubara summations over we in the ex- 
pressions for M11 an d M12 leads to 



c>qp-qh 



e 



n,ll 



M?i(q) = -g n 

m?m = 0rr 2 qh + 0^ p r2 qp . 



oqp-qp 

w n,ll i 



(25) 
(26) 



where the explicit form of the O n ,y functions is given in 
Appendix A, for the peruse of the reader. We choose 
to separate the contributions of the matrix elements 
Mij in terms of quasiparticle-quasiparticle (qp — qp) and 
quasiparticle-quasiholc (qp — qh) processes in order to 
isolate the channels that contribute to Landau damp- 
ing of the collective modes to be discussed in the next 
section. 



III. GAUSSIAN FLUCTUATIONS 



IV. COLLECTIVE MODES AT T = 



We now investigate Gaussian fluctuations in the pair- 
ing field $«(<?) about the the static saddle point A„.o- 
The Gaussian (quadratic) effective action can be written 
as 

S G = 5o(A„, ) + | £ A f (-g)M(g)A(g) (22) 
1 

where the fluctuation field is 

At(- ff ) = ( AI(g), Ai(- ? ), A 2 (-g), A*(g) ) (23) 
and M(g) is the fluctuation matrix given by 



M(g) = 



( M\i Mf 2 - 9 i2 

Mli M 22 -912 

-321 M| 2 M^ 

- ff 2i M 2 2 Mli 




where the matrix elements M-p^g) 
M? 2 (q)=M 2 \*(q) are given by 



= -g nn + (3- 1 Y,Wn( P + q/2)\ 2 MP + q)M-p) 
p 

M5 = /9- X X) r n(P + «/ 2 )«»(P + «)^(P)- 



The expressions i? n (p) = £ n (-p)/F(p) and ? n (p) = 
A n (k) /F(p) arc the matrix elements (G„.o)ii and 
(G n o)i2i respectively. Here, we use the definitions 
Up) = iwt ~ Cn(k), and F(p) = |£„(p)| 2 + |A„(k)| 2 . 
Notice that while Mi 2 (q) and M|i(g) are even under the 
transformations q — > — q and ive — > —ivi; MJj^g) and 
M" 1 (g) are even only under the transformation q — * — q, 
having no defined parity in we- In addition, notice that 
when gi2 = 321 = 0, and the fluctuation matrix becomes 



The collective modes are determined by the poles of the 
propagator matrix M _1 (g) for the pair fluctuation fields 
A(g), which describe the Gaussian deviations about the 
saddle point order parameter. The poles of M _1 (g) are 
determined by the condition detM = 0, and lead to a 
dispersion for the collective modes w — u>(q), when the 
usual analytic continuation ive — > w + i0 + is performed. 

We will focus here only at the zero temperature limit, 
but we will analyse phase and amplitude modes. At 
T = 0, only the qp — qp terms contribute, as the qp — qh 
terms vanish identically (see Appendix A). In this limit, 
we separate the diagonal matrix elements of M(g) into 
even and odd contributions with respect to w 



M?{ E {q) 



M?{°(q) = 



E 



T' 2 W + EE'}[E + E'] 
^ 2EE'[w 2 - {E + E') 2 Y { ' 

T' 2 [tt' + EE'}w 
2EE'[w 2 - {E + E') 2 }' 



(28) 



The off-diagonal term is even in w, and it reduces to 



m(o) = - 



T' 2 AA'\E + E'} 



^ 2EE'\w 2 

k 



(E + E*) 2 }' 



(29) 



We used in the previous expressions the following simpli- 
fied notation, the kinetic energies £ n (k) — > £, £„(k+q) — > 
the quasiparticle energies E n (k) — > E, E n (k + q) — > 
E'; the order parameters A„(k) — > A, A„(k + q) — > A'; 
and the symmetry factors T„(k) — > T, r„(k + q/2) — * T' . 

In order to obtain the collective mode spectrum, we 
express A n {q) = T n {q)e 1 ^ = (A„(g) + i0 n {q))/y/2 
where T n (q), <p n {q), X n (q) and n (q) are all real. Notice 
that the new fields A„(g) = r„(g) cos 0(g), and 9 n {q) = 
r„(g) sin 0(g) can be regarded essentially as the ampli- 
tude and phase fields respectively, when </>(g) is small. 
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This change of basis can be described by the following 
unitary transformation 



with W r , 



012 



x\X2- Again, the dispersion relation 




(30) 



Since we are considering the case where the saddle point 
order parameters Ai(k) and A2(k) are in phase the fluc- 
tuation matrix in the rotated basis then reads 



M 



where x n = M™{ 

rl,0 



Xi 


iz\ 





~.9l2 


—iz\ 


Hi 


-512 








-921 


2/2 


— iz 2 


-.921 





IZ2 


X2 



M 



12 J Vn 



(31) 



Mf 2 , and 



z n = Mi{ with the q dependence being implicit. The 
rotated matrix M is as complex as the un-rotated matrix 
M, and the real advantage of working in the rotated ba- 
sis involving fields A„ and 6 n is that the interpretation in 
terms of amplitude and phase fields is straight-forward. 
For instance, by inspection, notice that the phase fluctu- 
ation fields 9i and #2 corresponding to different conden- 
sates that are coupled by interband elements gi2 = 921- 
The same applies to amplitude fluctuation fields Ai and 
A2. When the interaction V12 = 0, the matrix element 
<?i2 vanishes, the matrix M becomes block diagonal and 
the two bands are decoupled. In this case the one-band 
results previously obtained are recovered for each inde- 
pendent band^l. 

Next, we focus on phase-phase and amplitude- 
amplitude collective modes. The easiest way to get the 
amplitude-amplitude collective modes is to integrate out 
the phase fields to obtain an amplitude-only effective ac- 
tion 

^ = f X>t( ? )Ma( g )A( g ), (32) 
1 

where \'(q) = (Xi(q), A2 (<?)). The amplitude-amplitude 
fluctuation matrix has the form 



M a = 



xi + V2z\lW a 

-512(1 - ZxZ 2 /W a ) 



-512(1 - ziz 2 /W a ) 

x2 + y\zl/W a 



(33) 



where W a = gf 2 — ViV2- The dispersion relation for the 
amplitude-amplitude collective modes is obtained from 
the condition detM a = 0. In this paper, however, we 
are mostly interested in the phase-phase collective modes. 
Thus, upon integration of the amplitude fields we obtain 
a phase-only effective action 



(34) 



where <%) = (e 1 (q),9 2 (q)). 
tion matrix has the form 

ivr - ( yi + x 2zf/w p 

p \ -912(1- z lZ2 /W p ) 



The phase-phase fluctua- 



-512(1 - ziz 2 /W p ) 
V2 +x 1 z^/W p 



(35) 



for the phase-phase collective modes is obtained from 
the condition dct M p = corresponding to poles of the 
phase-phase correlation matrix M p . In the next sec- 
tion, we discuss both analytical results for the phase- 
phase modes in BCS and BEC limits, as well as numerical 
results and the crossover regime. 

We note in passing that at finite temperatures, for con- 
tinuum s-wave systems, the qp — qp terms are well be- 
haved in the long- wavelength and low frequency limit, 
while the qp — qh terms do not in general allow for a 
simple expansion in the same limit due to the presence 
of Landau damping (see Appendix A). However, at zero 
temperature, the qp — qh terms vanish, and a well de- 
fined expansion is possible at low frequencies provided 
that the collective modes can not decay into the two- 
quasiparticle continuum. Thus, the collective mode dis- 
persion w(q) must satisfy the following condition 

iu(q) < min{£i(k)+£i(k + q),£;2(k) + £2(k + q)}. 

(36) 

To obtain the long wavelength dispersions for the collec- 
tive modes at T = 0, we expand the matrix elements of 
M (Eq. I|31l) ) in the amplitude-phase representation to 
second order in |q| and fourth order in w to get 



-*ln 

Pn 
Br, 



Qnfaf 

+ H, 



D n w 2 
R n w 2 



-E n w 2 \q_\ 2 
S n w 2 \q\ 2 - 



w 



- F n w\(37) 
T n w\(38) 
(39) 



where expansion coefficients arc given in Appendix B. As 
it will become clear in section [V] the |q| 4 order terms in 
the expansion are not necessary to calculate the collective 
mode frequencies w(q) accurately to order |q| 2 . 

The expressions of the coefficients found in Appendix 
B are valid for values of the interaction range parameter 
k n fi that satisfy the diluteness condition (k n fi/kF S> 1). 
However, in order to make analytical progress in the cal- 
culation of collective modes to be discussed next, we take 
the limit k n fi — > 00, since all the momentum integrals are 
convergent. This is in contrast to the situation encoun- 
tered with the order parameter equation, where we used 
k n ,o ~ 10 4 fci? in order to ensure convergence of our nu- 
merical calculations for |A„| and /1 as will be seen in the 
next section. 



V. ANALYTICAL AND NUMERICAL RESULTS 

In this section, we will focus only on the long- 
wavelength (small |q|) limit phase-phase modes deter- 
mined by the condition det M p = 0. We begin our dis- 
cussion with the trivial case when gi2 = (where V12 = 
but Vn and V22 can have any negative value) which corre- 
sponds to two uncoupled bands. In this case, the phase- 
phase fluctuation matrix M p = becomes block diag- 
onal, and we find two Goldstone modes satisfying the 
relation 



c^lql 2 , 



(40) 
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where the square of the speed of sound c„ in n band is 
given by 

c 2 = AnQn . (4i) 

B n + A n R n 

In the limit q — > 0, the corresponding eigenvectors are 
given by 

0% = O,w = O) = (0i,0a)oc(l,O) (42) 

for the first band, and 

0t(q = O,u; = O) = (0i,0 2 )oc(O,l) (43) 

for the second band, respectively. These results are iden- 
tical to the known results in the one-band caseSl. 

However, in the non-trivial case when gi 2 7^ 0, which 
corresponds to a finite Jopephson coupling Vyx between 
bands, we find two modes. In order to determine the 
collective mode spectra w(q) accurately to order |q| 2 it 
is sufficient to rewrite the determinant condition in the 
form 

w i (a 1 + a 2 |q| 2 ) + w 2 {a 3 + a 4 |q| 2 ) + a 5 |q| 2 = 0, (44) 

where a n are non-trivial and extremely complicated 
functions of the expansion coefficients A n ,B n , C n , D n , ... 
given in Appendix B. The exact dependence can be ob- 
tained via a symbolic manipulation program, but we will 
not quote these general results here or in appendices, as 
they are not particularly illuminating. Instead, we will 
present simple limits, where their behavior can be easily 
understood. 

In this non-trivial case of 1712 7^ 0, we find two collective 
modes. The first mode is the Goldstone mode satisfying 
the relation 

^ 2 -c 2 |q| 2 , (45) 

where the square of the speed of sound c is given by 

c 2 = ~a 5 /a 3 > 0. (46) 

In the limit q — > 0, the eigenvector of the Goldstone 
mode is given by 

0t(q = 0, w = 0) = (0i,0a) oc (|Ai|, |A 2 |) (47) 

which is valid for all values of the V\\ and V22 couplings. 
Notice that this mode is associated with the in-phase 
fluctuations of the phases of the order parameters around 
their saddle point values. 

In the particular case, where 312/ min{iVi, .N2} is the 
smallest expansion parameter (small 312 limit) we can 
perform a Taylor expansion of the a n coefficients around 
912 = 0, and obtain 

c 2 = t 1 t 2 (P 2 Q 1 + PiQ 2 )/(Piti + P 2 t 2 ), (48) 

as the square of speed c of the Goldstone mode. Here, 
we introduced a coefficient 

t n = (A n - P n )/[Bl + R n (A n - P n )} > 0, (49) 



which is positive definite since A n > P n and R n > 0. 
The precise meaning of the smallest expansion parame- 
ter 912/ min{iVi, N 2 } will be clear in sub-sections A and 
B, where we discuss analytically the weak and strong 
coupling limits. The eigenvector in the small 1712 limit 
is exactly the same as in the case of general (712, since 
0T(q = 0, w = 0) = (01,02) is independent of the value of 
312. (See Eq. gJJ) 

The second mode is a finite frequency mode satisfying 
the relation 

w 2 =w 2 +v 2 \q\ 2 1 (50) 
where the square of a finite frequency wq of the mode is 

wl = -a 3 / 'an > (51) 

and the square of the speed v of the mode is 

v 2 = a 3 a 2 /a\ - a^/ai + a 5 /a 3 . (52) 

The eigenvector for this mode has a complicated expres- 
sion 

f (q,w) = (01,0a) oc (gi 2 (W p - z x z 2 ), y\W p + x 2 z 2 ), 

for a general value of gi 2 . 

In the small gi 2 limit, the coefficients simplify to 

w 2 = Pih + P 2 t 2 > 0, (54) 
v 2 = Qih + Q 2 t 2 - c 2 > 0. (55) 

In the limit of w = wq, q — > 0, and small gi 2l the eigen- 
vector expression simplifies to 

6\q,w) = (01,0a) cx (|Aa|ti,-|Ai|t 2 ), (56) 

and it becomes transparent that this mode is associated 
with out-of-phase fluctuations of the phases of the order 
parameters around their saddle point values, since t n are 
positive definite. 

Before discussing the collective modes in the analyt- 
ically tractable weak and strong coupling limits, notice 
that a finite q Goldstone mode is only possible when 
c 2 > 0, and that a finite q finite frequency mode is only 
possible when Wq > 0. If these conditions are violated 
the modes are non-existent. Furthermore, caution should 
be exercised by recalling that the small w approximation 
used in the expansion of general fluctuation matrix (M) 
elements breaks down when the frequency w of any of the 
collective modes moves up into the continuum of two- 
quasiparticle states. Therefore, our results are strictly 
valid only for 

w < min{2Ei(k), 2£: 2 (k)}, (57) 
which corresponds to 

w<min{2|A 1 |,2|A 2 |} (58) 
in the weak coupling (BCS) limit, and to 

w « min{2VH 2 + |Ai| 2 , 2 VdMl + #o) 2 + |A 2 | 2 } 

(59) 

in the strong coupling (BEC) limit. With these condi- 
tions in mind, we discuss next the weak and strong cou- 
pling limits. 
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A. Weak Coupling Limit 

The s-wave weak coupling limit is characterized by 
the criteria fi > 0, fi > Eq, fi w sp ^> |Ai|, and 
fi — Eq ^> IA2I. Analytic calculations are particularly 
simple in this case since all integrals for the coefficients 
needed to calculate the collective mode dispersions are 
peaked near the Fermi surface (see Appendix B). In ad- 
dition, we make use of the nearly perfect particle-hole 
symmetry, which forces integrals to vanish, when their 
integrands are odd under the transformation £ — > — £. 
For instance, the coefficients that couple phase and am- 
plitude modes within a given band {B n and H n ) vanish. 
Thus, in this case, there is no mixing between phase and 
amplitude fields within n th band, as can be seen by in- 
spection of the fluctuation matrix M. 

We would like to focus on the phase-phase collective 
modes, as they correspond to the low energy part of the 
collective mode spectrum. Notice that, all expansion co- 
efficients appearing in the phase-phase fluctuation matrix 
M p are analytically tractable. The expansion of the ma- 
trix elements to order |q| 4 and w 4 is performed under the 
condition (w, |q| 2 /2m„) < min{2|Ai|, 2|A 2 |}. To evalu- 
ate x n for each band n, we need A\ = gi 2 |A 2 |/|Ai| + N\, 
A 2 — (72i|Ai|/|A 2 | + N 2 , which are the coefficients of 
the (q = 0,w = 0) term; C n = c 2 >w N n /12\ A„| 2 , which 
are the coefficients of |q| 2 ; D n — A„/12|A„| 2 , which are 
the coefficients of w 2 ; and F n — — iV n /120| A„| 4 , which 
are the coefficients of w 4 . To evaluate y n , we need Pi = 
012 1 A 2 1 /| Ai| and P 2 = g 2 i|Ai|/|A 2 |, which are the coef- 



A simple evaluation of det M p = leads to a Goldstonc 



ficients of the (q = 0, w = 0) term; Q r< 
which are the coefficients of |q| 2 ; i?„ 
are the coefficients of w 2 ; and T n = 
are the coefficients of w 4 . Here, 



V n ,F 



/Vd„ 



,N n /A\A % 
= A n /4|A„| 2 , which 
-7V„/24|A„| 4 which 



(60) 



is the velocity of the sound mode in the one-band casei£, 
d n is the dimension, v n ,F is the Fermi velocity, and N n 
is the density of states at the Fermi energy per spin 
(N n = m n L 2 /2n in 2D and N n = m n L 3 k n _F /27r 2 in 
3D) in n th band. While we use N n as the density of 
states per spin at the Fermi energy, in Ref^ the density 
of states used includes both spins. The off-diagonal ele- 
ments are (M p )i 2 = (M p ) 2 i = -312, since B x = B 2 = 
and Hi — H 2 — 0, as discussed above. Notice that the 
expressions above are valid for both 2D and 3D bands. 

In order to bring our results in contact with Leggett's^, 
and Sharapov et. al£, we make use of the order param- 
eter saddle point Eq. (JT7I) at T = 

lAiKl + VuFi) = -\A 2 \V 2 iF 2 , (61) 
|A 2 |(1 + F 22 F 2 ) - -|Ai|Vi2*i, (62) 

and consider the small g\ 2 limit where 

3i 2 /min{JVi, N 2 } < min{|Ai|, |A 2 |}/max{|Ai|, |A 2 |}. 

(63) 



mode 



|q| 2 , characterized by the speed of sound 



N 2 c 2 



Ni+N 2 

and a finite frequency (Leggett) mode w 2 
characterized by 



(64) 



^|q| 2 , 



II u 



m+N 2 |8Vi 2 ||Ai||A 2 | 
2NiN 2 VnV 22 -V 2 2 
Nicj w + N 2 cj w 
N!+N 2 



(65) 
(66) 



where wq is a finite frequency, and v is the speed of prop- 
agation of the mode. Here we reintroduced all the cou- 
pling constants of the original Hamiltonian in Eq. ifljl. 
These results are valid only in the weak-coupling limit, 
with all V nrn < 0, and detV > 0. Notice that if V\ 2 = 
the Leggett mode does not exist as the two bands are 
uncoupled. Furthermore, the trivial limit of one-band 
(say only band 1 exists) is directly recovered by taking 
|A 2 | = 0, N 2 = and c 2jtu = which leads to c 2 — c\ wl 
w 2 = 0, and v 2 = 0. 

It is also very illustrative to analyse the eigenvectors 
associated with the two-solutions. For Goldstone's mode, 
in the limit of q — > 0; for any g\ 2 , it is easy to see that 



0t( q = O, W = O) = (0 1 ,^)cx(|A 1 |,|A 2 |), 



(67) 



corresponding to an in-phase mode. In the degenerate 
case, where |Ai| = |A 2 |, and N\ = N 2 , the eigenvector 
6^(q = 0, W = 0) oc (1, 1) and the phase fields are per- 
fectly in phase. The eigenvector for the finite frequency 
mode is T (q, iu) = {6-1,62) oc (gi 2 ,yi). In the particular 
limit where g\ 2 is small (Leggett mode), this leads to an 
eigenvector 



6> f (q = 0,w 



w ) 



} 6 2 )<x(N 2 \A 1 \,-N 1 \A 2 \), 



(68) 

which corresponds to an out-of-phase mode. In the de- 
generate case, this simplifies further to T (q = O, iu = 
0) = (1, —1) becoming a perfectly out-of-phase mode. 

Now, we would like to turn our attention to the analysis 
of the phase-phase modes in the strong coupling limit, 
which is also analytically tractable. 



B. Strong Coupling Limit 

The s-wave strong coupling limit is characterized by 
the criteria fx < 0, and |/z| » |Ai| and |/z| + E > |A 2 |. 
The situation encountered here is very different from 
the weak coupling limit, because one can no longer in- 
voke particle-hole symmetry to simplify the calculation 
of many of the coefficients appearing in the fluctuation 
matrix M (see Appendix B). In particular, the coeffi- 
cients B n 7^ indicate that the amplitude and phase 
fields within an individual band n arc mixed. 
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We will concentrate here only on phase-phase modes 
characterized by the fluctuation matrix M p . The ex- 
pansion of the matrix elements to order |q| 2 and w A is 
performed under the condition (w, |q| 2 /2m„) <C 2|/i|. 
All the coefficients (x n , y n , z n , W p ) appearing in M p ma- 
trix Eq. (|3l>jl are evaluated to order (|Ai|/|/z|) 2 and 
(|A 2 |/(|/i| + Eq)) 2 in the strong coupling limit for both 
two and three dimensional systems. To evaluate x n for 
each band n, we need A\ — gri 2 1 A 2 1/ 1 Ai | + ki|Ai| 2 /2|^|, 
A 2 = 52i|Ai|/|A 2 | + k 2 |A 2 | 2 /2(|^| + Eq), which are the 
coefficients of the (q = 0, w = 0) term; C n — K n /4m n , 
which are the coefficients of |q| 2 ; D\ — Ki/8|/i|, D 2 = 
K 2 /8(|/i| + Eq), which are the coefficients of w 2 , and 
Fi = -Ki7i/512|/i| 3 , F 2 = -k 27 i/512(|//| + E ) 3 , which 
are the coefficients of w A . To evaluate y n , we need 
Pi = <7i 2 |A 2 |/|Ai| and P 2 = t? 2 i|Ai|/|A 2 |, which are the 
coefficients of the (q = 0,w = 0) term; Q n = K n /4m n , 
which are the coefficients of |q| 2 ; R\ — Ki/8|/i|, R 2 = 
K 2/8(|/i| + Eq), which are the coefficients of w 2 , and 
Ti = - Kl7 /5f2|^| 3 , T 2 = -k 27 /512(|m| + E ) 3 , which 
are the coefficients of w 4 . To evaluate z n for each band 
n, we need B n = k„, which are the coefficients of w, and 
Hi = K17/96H 2 , H 2 = K 2 7/96(|/i| + Eq) 2 which are the 
coefficients of w 3 . The computation of W p — g\ 2 — X\x 2 , 
just relies on the knowledge of x n already obtained above. 
While the expressions above are valid for any value of /j, 
and Eq in the 2D bands, they are only rigorously valid 
for ^> Eq in the 3D case through the use of the di- 
mensionally dependent variables K n , and 7 and 7 which 
in the 3D case become 



«i = TrNi/8-\/\Jj^, 

k 2 = ttN 2 /8^(\ij i \ + E )s F , 



(69) 
(70) 



and 7 = 5, 7 = 3. To write our expressions in compact 
notation and to recover the degenerate limit (Eq — > 0), 
we keep Eq in all 3D coefficients. All expressions for the 
3D case can be converted into the 2D case through the 
following procedure. In all coefficients (x n ,y n , z n ,W p ), 
the dependence in + Eq) is transformed into 

|/x|/2, + Eq)/2, and the variables K n , and 7 and 7 
need to be redefined as 



K2 



Ni/8\n\, 
N 2 /8{\fx\+Eo), 



(71) 
(72) 



and 7 = 2, 7 = 2. 

In the limit q — * 0, the condition detM p = leads 
again to two modes. General expressions for the col- 
lective modes arc highly non-trivial and, therefore, we 
consider analytically only the asymptotic small g 12 limit 
defined as 



W + Eq)} 2 « 1. 

(73) 

w In this case, the Goldstone mode corresponds to w 2 — 

„2|„|2 



ffl2 /min{Ar 1 ,7V 2 }< [|min{|A 1 |,|A 2 

w In this cast 
c 2 |q| 2 , where 



Kl\^\cl s + K 2 (\fl\ + E )4 S 

Ki\(jl\ + k 2 (\li\ + Eq) 



(74) 



is the square of the speed of sound. The finite frequency 
mode (which is the extension of Leggett's mode in the 
weak coupling limit) corresponds to w 2 — w 2 |q|' 
where 



^0; 



«iM + « 2 (H + £o) [VigAiAaj 
2 Kl \fi\ K2 (\n\ + Eq) VnV 22 -V 2 
Ki\y\cj s + k 2 {\/i\ + E )cl s 

Ki\/J,\ + K 2 (\H\ + Eq) ' 



(75) 



(76) 



are the finite frequency, and speed of propagation of the 
mode, respectively. Here, the quantitites 



c M = |Ai|/V8mi|/i|, 
c 2 . s 



\A 2 \/y/8m 2 \(fi\ +Eq) 



(77) 
(78) 



are the velocities of the sound mode in the one-band 
case2i. These results are valid only in the strong-coupling 
limit, with all V nm < 0, and detV > 0. Notice that if 
Vi 2 = the finite frequency mode does not exist as the 
two bands are uncoupled. Furthermore, the trivial limit 
of one band (say only band I exists) is directly recovered 
by taking k 2 — and C2, s — which leads to c 2 = c\ s , 
w 2 = 0, and v 2 = 0. 

The eigenvectors associated with these solutions are as 
follows. For Goldstone's mode, in the limit of q — ► and 
for any value of g\ 2 , 



0t( q=O;U , = o) = (6 U 6 2 ) oc (|Ai|,|A 2 |), 



(79) 



corresponding to an in-phase mode. For the finite fre- 
quency mode (in the small g\ 2 limit), the eigenvector 
becomes 

e^(q = 0,w = wo) cx (/safl/il + ^WIAil.-Sil/illAal), 

which corresponds to an out-of-phase mode, as k„ > 0. 
In the degenerate case (k\ = k 2 , |Ai| = |A 2 | and Eq — > 
0), this simplifies further to 6^(q = 0, w = 0) = (I, — I) 
becoming a perfectly out-of-phase mode, similar to the 
weak coupling case. 

Next, we would like to turn our attention to the anal- 
ysis of the phase-phase modes in the crossover region 
which is not analytically tractable and requires numeri- 
cal calculations. 



C. Numerical Results and Crossover Region 

Thus far, we have focused on the analytically tractable 
limits corresponding to weak and strong couplings. In 
order to gain further insight into the behavior of the 
phase-phase collective excitations at T — 0, we present 
numerical results of the evolution from weak to strong 
coupling for the Goldstone and finite frequency modes. 
We limit ourselves to numerical calculations of the fully 
degenerate (identical) bands case, from which the known 
one-band results^ for the Goldstone mode can be easily 
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recovered in the BCS, BEC, and crossover regimes for an 
s-wave supcrfiuid. While the limit of degenerate bands 
is probably harder to find in nature, it provides us with 
qualitative and quantitative understanding of the evolu- 
tion of the finite frequency collective modes from weak to 
strong coupling. At the same time the degenerate prob- 
lem is easier to solve numerically and serves as a test 
model to our analytical results. Thus, we postpone a 
detailed numerical calculation of non-denegerate bands 
case for a future publication, where the more complex 
numerical problem will be attacked. 

In this particular case (identical bands), the band offset 
is Eq — > 0, the density of states at the Fermi energy are 
N% = N 2 = N; the Fermion masses are mi — m 2 = m; 
the Fermi velocities are v\^f = v 2 ,f = vf', and the 
intra-band interactions are V\\ = V 22 — V; while the 
order parameter amplitudes are |Ai| = IA2I = |A|). 
Furthermore, the low-momentum and low-energy expan- 
sion coefficients become identical, i.e., A\ — A 2 = A: 
B 1 = B 2 = B- D l = D 2 = D; H x = H 2 = H; 
Ri = R-2 = R; and T\ = T 2 = T, and we define 
A = A - 5 i2 = J2k A(k) 2 /2£(k) 3 . For any value of 
(712 in the degenerate bands limit, we obtain 



QA 



B 2 +AR 



>0, 



(81) 



for the squares of the speeds c (Goldstone) and v (finite 
frequency). In addition, we obtain 



m 9 



1.2) 



(82) 



for the square of the finite frequency which is also valid 
for any value of g± 2 . The numerator V 2 {g\ 2 ) = 2{B 2 + 
AR)(gi 2 A + 2g\ 2 ) is positive definite and the denomina- 
tor P 2 (g 12 ) = (B 2 +AR) 2 +g 12 [R(B 2 + AR)+A(2TA~ 
4BH) + 8ARD + 6B 2 D] + g 2 l2 [4TA + 8{RD - BH) + 
5B 2 D/A] must be also positive definite to guarantee that 
Wq > 0. Notice that these functions are not strictly sec- 
ond order polynomials in gi 2 , since all coefficients de- 
pend implicitly on 1712. In the small gi 2 limit, expressions 
for c 2 , v 2 and w 2 simplify to Eqs. (|55jl and 

with the degenerate case coefficients. Furthermore, since 
'^'2(512) > is positive definite, Eq. (|82[1 is only valid 
strictly for V 2 {gi 2 ) > 0. 

In our numerical calculations of |A„| and fj, (via the 
saddle point Eqs. (|17|l and \21$ ). we choose a momentum 
cut-off of value fei = ^2,0 — ^0 ~ 10 i kp to ensure con- 
vergenge of all k-space integrations. The fermion density 
can be expressed in dimensionless units as n s = n/n 
where n max = k F 2 Jlax /2Tr in 2D, and n max = k F max 
in 3D. The maximal value of the Fermi momentum kp max 
(that fixes the maximal density n ma x) is chosen by fixing 
the ratio kpmax/ko — 10~ 4 , which easily satisfies the di- 
luteness conditions (ko/kp max ) 3 ^ 1 ( or n maxRo <C 1) 
in 3D, and {k /k Fmax ) 2 > 1 (or n max R% < 1) in 
2D. See Sec. [H] to recall how the interactions depend 



max , 

k F 3 ma x/^ 2 



on ko = 27r/i?o, where Rq plays the role of the inter- 
action range in real space. For any kp /kFmax < 1 
all conditions are satisfied, thus we work at fixed den- 
sity n s = 1/27T w 0.159, corresponding to kF/kp max = 
1/V2n = 0.398 in 2D, and n s = 3/(8tt) w 0.119, corre- 
sponding to k F /k Fmax = (3/871-) 1 / 3 w 0.492 in 3D. 

We also confine ourselves to the asymptotic small g\ 2 
limit, which means g\ 2 /N <C 1 in weak coupling, and 
9n/N <C [|A|/|/i|] 2 in strong coupling limits. Therefore, 
in 3D, we choose V\ 2 = 10 _7 V (since V\\ — V 22 = V) 
which leads to gi 2 /N ~ 10 -4 , for a 3D Fermion density 
of n s — 0.119. In the 2D case, we choose V12 = 10 — 5 V^, 
which leads to gi 2 /N ~ 10 -4 , for a 2D Fermion density 
of n s = 0.159. This particular choice satisfies the small 
gi 2 condition for the range of couplings shown in Figs. 2 
and 4. 

We solve the saddle point equations for the order pa- 
rameter Eq. (|17|) together with the number equation 
Eq. H21fl self-consistently for fixed densities. The order 
parameter amplitude |A| and chemical potential fi in 2D 
(3D) are presented in Figs. 2a and 2b (Figs. 4a and 4b) 
as a function of the dimensionless intra-band interaction 
parameter V r — NV/ir. Notice that the system crosses 
over from the BCS (fj, > 0) to BEC (/i < 0) regimes at 
V r = 1.33 x 10~ 2 in 2D, and at V r = 1.59 x 10~ 4 in 3D, 
where the chemical potential // crosses zero (the bottom 
of the degenerate bands). 

For the 2D (3D) case, we show in Figs. 3a and 3b 
(Figs. 5a and 5b), numerical plots of the sound velocity 
c, normalized by the Fermi velocity vf, and of the ratio 
of the finite frequency wq with respect to the minimum 
quasi-particle excitation energy min{2i?(k)} as a func- 
tion of intra-band couplings V r . In addition, notice the 
very good agreement between the numerical results and 
the analytical approximations in their respective (BCS 
or BEC) limits. 

The analytical value for the weak coupling sound ve- 
locity follows from Eq. (|64|l for the non-degenerate case 
with Ni — N 2 and c\ )W = c 2lW , which leads to c 2 = vp/2 
for 2D and c 2 — v F /3 for 3D bands. Similarly, the analyt- 
ical value for the strong coupling sound velocity follows 
from Eq. (|74H with K\ = k 2 and ci iS = c 2iS , which leads 
to c 2 = |A| 2 /4m|^| for 2D and c 2 = |A| 2 /8m|^| for 3D 
bands. Thus, we recover the Goldstone mode in both 
BCS and BEC limits as in the case of the presence of 
only one band 1 ^!. The numerical values (solid circles 
in Figs. 3a and 5a) for the sound velocity as a func- 
tion of the dimensionless coupling V r are calculated from 
Eq. (18 1 11 . Notice the very good agreement with the ana- 
lytical results in weak and strong coupling (dotted lines 
in the same figures). As a further consistency check, no- 
tice that this agreement is very reasonable since Eq. I|81|) 
is identical to the expression for the sound velocity given 
in Refi^£ for the one-band model, with the correspon- 
dence that our coefficient A has the same expression as 
the coefficient A defined in their paper. 

Notice in Fig. 3a that the sound velocity c is essentially 
a constant for all couplings V r m the ko/kFmax ~ * 00 
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FIG. 2: Plots of a) chemical potential n scaled by ef and b) order parameter |A| scaled by ef versus dimensionless coupling 
V r (see text for definition) for two-dimensional degenerate bands. The region where fj, changes sign is shown in the inset. Note 
that fj, = when V r = 0.0133. 
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FIG. 3: Numerical calculation (solid circles) of a) square of sound velocity c 2 scaled by v F , and of b) the square of finite 
frequency w 2 scaled by the two-quasiparticle threshold min{2_B(k)} versus dimensionless coupling V r for two-dimensional 
degenerate bands. The dotted lines represent analytical results for the weak and strong coupling limits. Notice that, we scaled 
Wo by 4|A| 2 for fj, > 0, and by 4(|^| 2 + |A| 2 ) for /x < 0. 



(ko/kprnax = W ) limit. For smaller values of ko/kF max 
the sound velocity decreases as a function of V r (not 
shown in figure). We will not discuss here the depen- 
dence of c on the ratio ko/kF ma x since we are mostly 
concerned with checking the consistency of our calcula- 
tions with the analytically tractable limits. In this case, 
the sound velocity c is not a good indicator that the BCS- 
BEC crossover is occuring in 2D. However, in the 3D case 
(see Fig. 5a), the speed of sound changes very fast in the 
neighboorhood of fi = 0, thus manifesting itself as an in- 
dicator of the crossover regime. Therefore, measurements 
of Goldstone mode frequency can offer an indication of 
the BCS-BEC crossover possibly only in 3D two-band 
superfluids. 

With the confidence of recovering the sound mode re- 
sults for a one-band model from a two-band model with 
identical bands, we proceed with the discussion of the fi- 
nite frequency mode, which is shown in Figs. 3b and 5b 
for the 2D and 3D cases, respectively. 

The analytical value of Wq for the weak coupling fi- 
nite frequency mode follows from Eq. I|b5|) for the non- 
degenerate case with Vn = V22 = V, Ni = N 2 = N, and 



|Ai| = |A 2 | = |A|, which leads to w 2 = 8gi 2 |A| 2 /iV 
for both 2D and 3D bands. Similarly, the analyt- 
ical value for strong coupling follows from Eq. J7SJ), 
which leads to w 2 — 8gi2|A| 2 /iV for 2D, and w% = 
8gi2|A| 2 \/£i?/|//|/7riV for 3D bands, respectively. Nu- 
merical values (solid circles in Figs. 3b and 5b) for the 
finite frequency wq as a function of the dimensionless cou- 
pling V r are calculated from Eq. I|82|) . Notice the very 
good agreement between the numerical results and their 
analytical counterparts in both weak and strong coupling 
limits (dotted lines in the same figures). It is important 
to notice that the scales for weak and strong coupling- 
used in the finite frequency plots are not the same. We 
scaled Wq by min{2£ , (k)}, which corresponds to the two- 
quasiparticle excitation threshold, thus Wq is scaled by 
4|A| 2 since min{2£:(k)} = 2|A| for fx > 0, wl is scaled 
by 4(|/x| 2 + |A| 2 ) since min{2£(k)} = 2 V ^| 2 T|A| 2 
for fi < 0. This choice is natural, because it indi- 
cates that the finite frequency wo always lies below the 
two-quasiparticle excitation threshold for the parameters 
used, meaning that the collective is undamped. 

Notice in Figs. 3b and 5b {ko/kF max —> 00, 
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FIG. 4: Plots of a) chemical potential n scaled by ef and b) order parameter |A| scaled by ef versus dimensionless coupling 
V r (see text for definition) for three-dimensional degenerate bands. The region where [i changes sign is shown in the inset. 
Note that \i = when V r = 1.592 x 10" 4 . 
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FIG. 5: Numerical calculation (solid circles) of a) square of sound velocity c 2 scaled by v F , and of b) the square of finite 
frequency Wo scaled by the two-quasiparticle threshold min{2_B(k)} versus dimensionless coupling V r for three-dimensional 
degenerate bands. The dotted lines represent analytical results for the weak and strong coupling limits. Notice that, we scaled 
Wo by 4|A| 2 for fj, > 0, and by 4(|/i| 2 + |A| 2 ) for fj, < 0. 



ko/kFmax = 10 4 ) that the finite frequency wq changes 
qualitatively near the coupling V r where [i changes sign 
for both 2D and 3D cases. We will not discuss here 
the dependence of wq on the ratio kQ/kp max since we 
are mostly concerned with checking the consistency of 
our calculations with the analytically tractable limits. 
However, it is important to emphasize that the finite 
frequency mode is a good indicator that the BCS-BEC 
crossover is occuring in both 2D and 3D. Thus, measure- 
ments of the finite frequency wq can reveal BCS-BEC 
crossover behavior in two-band superfluids. 



VI. CONCLUSIONS 

We studied the evolution of low energy collective exci- 
tations from weak (BCS) to strong (BEC) coupling limits 
in two-band s-wave superfluids at T = for all intra-band 
coupling strengths with ranges satisfying the diluteness 
condition. We assumed that the two bands were coupled 
via an inter-band Josephson interaction. We focused on 
the phase-phase collective modes and showed that there 



can be two undamped phase-phase modes in the evolu- 
tion from weak to strong coupling. In the weak coupling 
limit, we recovered Leggett's results corresponding to an 
in-phase mode (Goldstone) mode and an out-of-phase 
mode (Leggett's mode) in the appropriate asymptotic 
limits. Furthermore, we generalized Leggett's weak cou- 
pling results to include the BCS-BEC crossover and the 
strong coupling regime. In addition, we presented analyt- 
ical results in the strong coupling limit, in the asymptotic 
limit of small gi2 corresponding to weak Josephson cou- 
pling between the bands. All the analytical results were 
presented for the cases of two-dimensional and three- 
dimensional bands, and the for cases of non-degenerate 
and degenerate bands. 

On the numerical side, we analysed fully the limit 
of degenerate bands, from which the one-band results 27 
can be easily recovered in the BCS, BEC, and crossover 
regimes for a 3D s-wave superfluid. The limit of degen- 
erate bands, although less likely to be found in nature, 
provides us a good basis for the more challenging numer- 
ical work for the non-degenerate case, which will be per- 
formed in the future, as they require the self-consistent 
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solutions of three simultaneous non-linear integral equa- 
tions in order to determine the order parameters |Ai|, 
IA2I, and the chemical potential fi. 

We have also briefly described in this manuscript the 
procedure to compute the amplitude-amplitude collective 
modes, although we have not discussed in detail their ex- 
plicit form, we would like to mention that amplitude- 
amplitude collective modes are higher energy modes 
in comparison to lowest energy (phase-phase) collective 
modes discussed here. The phase-phase and amplitude- 
amplitude collective modes are potentially important to 
the analysis of multi-component ultracold neutral Fermi 
gases, where collective mode frequencies can be measured 
spectroscopicall y 24 ! 25 . 

The results for collective modes described here (neu- 
tral) is not strictly valid in standard (charged) condensed 
matter systems. In particular, in the BCS limit, the effect 
of the Coulomb interaction is to plasmonize the the Gold- 
stone mode. Furthermore, the Coulomb interaction mod- 
ifies the velocity of the Leggett mode, but does change 
its finite frequency offset^. There is some experimen- 
tal evidence that the Leggett mode in the BCS limit has 
been observed^ in MgB2, but there is presently no ex- 
perimental charged two-band condensed matter system 
where the BEC limit can be reached via the tunning of 
an experimental parameter. The interesting extension to 
the BEC limit and the crossover region in the charged 
case is currently underway and will be published else- 
where. 

To conclude, the main contribution of our manuscript 
is to study the evolution of the Goldstone (in-phase) 
and the finite frequency (out-of-phase) collective modes 
from weak (BCS) to strong (BEC) couplings in neutral 
two-band superfluids. Our results are potentially rele- 
vant to multi-component ultra-cold Fermi atoms, and can 
be used, in principle, to test the validity of the BCS- 
BEC evolution based on extensionsi£ii& of Eagles—, 
Leggett'si^, and Nozieres and Schmitt-Rink'si& sugges- 
tions. 
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In the previous expressions the indices qp — qh 
and qp — qp classify quasiparticle-quasihole and 
quasiparticle-quasiparticle terms, respectively. Fur- 
thermore, we used the following simplified notation: 
the kinetic energies £ n (k) — > £, £ n (k + q) — > the 
quasiparticle energies E n (k) — > E, E n (ln + q) — > E'; the 
order parameters A„(k) — > A, A n (k + q) — » A'; the 
symmetry factors r„(k) — > T, r„(k + q/2) — * T'; and the 
Fermi functions f n (E n (k)) -> /, f n (E n {k + q)) -> /'. 
We also made use of the definition of the first coherence 
factor 

|h| 2 = ~(1 + |), (Al) 
[\u'\ 2 = (1 + £' /E')/2], the second coherence factor 

M 2 = ^(l-f) (A2) 



[\vf 
them 



(1 — £'/E')/2], and the phase relation between 



A 

2E 



(A3) 



[u'*v' = A'/2E']. Finally, we used the notation X = 
if - /') r ' 2 and Y = (1 — f — f')T' 2 to indicate the com- 
binations of Fermi functions (/, /') and symmetry coef- 
ficients (r, r') appearing in the quasiparticle-quasihole, 
and quasiparticle-quasiparticle terms, respectively. 
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APPENDIX A: MATRIX ELEMENTS 

In the evaluation of the elements of the fluctuation 
matrix M(g) appearing in section IIIII and defined in 
Eqs. ffiijl. we need to calculate the functions 6 I | P n qh , 
©^Li qP > e n P r 2 ql \ and 6^r 2 qP given by 



k 



MVI 2 



iv + E- E' 



iv 



APPENDIX B: EXPANSION COEFFICIENTS 

From the rotated fluctuation matrix M expressed in 
the amplitude-phase basis as defined in section llVl we can 
obtain the expansion coefficients necessary to calculate 
the collective modes at T = 0. In the long- wavelength 
(q — > 0), and low frequency limit (w — > 0) the matrix M 
defined in Eq. (|31f) is fully determined by the knowledge 
of the matrix elements x n , y n and z n . Note that, this 
expansion requires 



7, \q\ 2 /m n « min{^ 1 (k),S 2 (k)} 



(Bl) 



and these coefficients are valid for all couplings for s-wave 
pairing. 

In all the expressions below we use the following sim- 
plifying notation f = dT/dk, £ = d£/dk, f = d 2 T/dk 2 , 
£ = d 2 £/dk 2 , and k = |k|. 
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The coefficients necessary to obtain the matrix element 



An 



e 



V^r 2 



k 



(B2) 



corresponding to the (q — 0, w — 0) term, 

c « = E^{^ 2 - 3A2 ) r2 + ^ 2 ( 2A2 -^ 2 ) r r 



2« 



- I OA 2 ( 1-^2 
+ A 2 ( I. - 10-^ 



T 2 cos 2 a 



£ 2 (£ 2 -7A 2 )-5A 4 x _2 



FT cos 2 a (B3) 



r 2 cos 2 



corresponding to the |q| 2 term with a being the angle 
between k and q; 



8# 5 



corresponding to the w 2 term, and 



« 2 



f - V g r 



(B4) 



(B5) 



corresponding to the w 4 term. 

The coefficients necessary to obtain the matrix element 
Vn are 



Pn 



V— r 2 

Z-2F ' 



corresponding to the (q — 0,w — 0) term, 

Q n = J2 8^5 {^ 2r2 - ^ 2 - 3A2 ) r2 - t 2 E 2 rf 

k 

[ti(2E 2 - 6A 2 )rr - £ 2 (£ 2 - 2A 2 )F 2 ] cos 2 aj(B7) 



+ 



corresponding to the |q| 2 term with a being the angle 
between k and q; 



r 2 

8£ 3 ' 



corresponding to the w 2 term, and 



(B8) 



T n = -Y — 



32E 5 



(B9) 



corresponding to the w 4 term. 

The coefficient necessary to obtain the matrix element 



z„ is 



€ -p2 
4S 3 ' 



corresponding to the w term, and 



i p2 



16£ 5 



(BIO) 



(Bll) 



(B6) corresponding to the w 3 term. 
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